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ABSTRACT 

We estimate the variance £ 2 ; the skewness ^ 3 and the kurtosis £ 4 in the distribution of 
density fluctuations in a complete sample from the APM Cluster Redshift Survey with 
339 clusters and a mean depth V ~ 250 Mpc* We are able to measure the statistics 
of fluctuations in spheres of radius R ~ 5 — 80 Mpc, with reasonable errorbars. The 

statistics in the cluster distribution follow the hierarchical pattern £j = Sj ^ 2 with 
Sj roughly constant, S3 ~ 2 and S4 ~ 8. We analyse the distribution of clusters taken 
from N-body simulations of different dark matter models. The results are compared 
with an alternative method of simulating clusters which uses the truncated Zel'dovich 
approximation. We argue that this alternative method is not reliable enough for making 
quantitative predictions of £. The N-body simulation results follow similar hierarchical 
relations to the observations, with Sj almost unaffected by redshift distortions from 
peculiar motions. The standard SI = 1 Cold Dark Matter (CDM) model is inconsistent 
with either the second, third or fourth order statistics at all scales. However both_a 
hybrid Mixed Dark Matter model and a low density CDM variant agree with the £j 
observations. 



Key words: Large-scale 
numerical, statistical. 



structure of the 



universe 



galaxies: clustering, methods: 



INTRODUCTION 



Rich clusters of galaxies can be used to trace the structure 
^5 °f the Universe on very large scales. Although the identi- 



C3 



fication of clusters is more difficult than that of galaxies, 
clusters have several important advantages as tracers, both 
in observations and in simulations. Since clusters represent 
high densities in the galaxy distribution, it is possible to 
use them to map out the distant Universe in a systematic 
and economical way. As clusters are associated with high 
peaks in the mass distribution, cluster correlations evolve 
very slowly with time, and so the clustering statistics deter- 
mined from numerical simulations are relatively insensitive 
to the epoch at which the simulations are analysed (Croft & 
Efstathiou 1994a). An additional advantage is that clusters 
are strongly correlated so that even a small sample can be 
used to estimate the correlation functions. 

The 2-point correlation function for galaxy clusters was 
first analysed in Abell's (1958) catalogue (e.g. Hauser & 
Peebles 1973, Bahcall & Soneira 1983, Klypin & Kopylov 
1983). Recent determinations based on automated cluster 
catalogues (Dalton et al. 1992, Nichol et al. 1992) and on flux 



Throughout the paper we use Hq = 100/i km/s/Mpc. 



limited samples of X-ray clusters (Romer et al. 1994) give 
mutually consistent results. The 2-point correlation of the 
APM (Automatic Plate Measuring machine) Cluster Red- 
shift Survey has already been used to constrain dark matter 
models (Dalton et al. 1994b). 

Previous analysis of 3-point and 4-point statistics, {3 
and {4, of cluster samples have been based on Abell's cat- 
alogue and extensions by Abell, Corwin & Olowin (1989) 
(e.g. Jing & Zhang 1989, Toth, Hollosi & Szalay 1989, Jing 
& Valdarnini 1991, Plionis & Valdarnini 1995, Cappi & Mau- 
rogordato 1995). These samples seem to reproduce the hi- 
erarchical relation { j ~ f 2 -1 found in the galaxy distribu- 
tion (e.g., Groth & Peebles 1977; Fry & Peebles 1978; Sharp 
et al. 1984; Szapudi, Szalay & Boschan 1992; Meiksin, Sza- 
pudi & Szalay 1992; Bouchet et al. 1993; Gaztanaga 1992, 
1994). These observations could provide important insights 
into the underlying matter distribution, and the formation of 
clusters. However, there is evidence that the clustering mea- 
sured from Abell's catalogue is affected by inhomogeneities 
in the selection of clusters (e.g. Efstathiou et al. 1992) and 
it is not clear how this could affect the above mentioned 
estimates of {3 and {4. 



Here we estimate the variance f 2 



the skewness f 3 



and 



the kurtosis f 4 in the distribution of density fluctuations 
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in a sample from the APM Cluster Redshift Survey. The 
selection of clusters in the APM Cluster Catalogue is based 
on a computer algorithm applied to digitised photographic 
plates and so provides a more uniform sample of clusters 
than Abell's. We compare the results with previous analysis 
and with different cluster simulations. 

In next section, we review the APM cluster data and 
N-body models. We present the method of estimation and 
the results in Section 3. Sections 4 and 5 are devoted to 
the discussion and conclusions. In an appendix we consider 
the properties of simulations generated using the truncated 
Zel'dovich approximation. 



sen to be 0.5 /t _1 Mpc and the mean intercluster separation 
30 /t _1 Mpc (Croft & Efstathiou 1994a). 

We have also tested an alternative scheme for mak- 
ing theoretical predictions about cluster clustering. This 
scheme, advocated by Plionis et al. (1995) makes use of the 
Zel'dovich (1970) approximation to dynamically evolve the 
density field before clusters are selected. In order to assess 
the validity of results derived using this method we have 
examined clustering in the Zel'dovich approximation and 
carried out a comparison with N-body clusters. These tests 
are detailed in Appendix A. 



2 CLUSTER SAMPLES 

2.1 The APM Cluster Survey 

We use sample B from Dalton et al. (1994b) which contains 
364 clusters with mean redshift 0.086 and space-density 
3.4 x 10 -5 /t 3 Mpc -3 . The clusters were selected from the 
APM Galaxy Survey (Maddox et al. 1990b, c) using an au- 
tomatic selection algorithm described in detail by Dalton 
et al. (in preparation). We have cut the sample to redshifts 
larger than lOOOOkm/s to avoid uncertainties in the selec- 
tion function at low redshift, and imposed an upper limit of 
40000km/s. There are only 25 clusters outside this range, so 
that the total number in the selected region is N = 339. As 
the redshifts are obtained from only a few galaxies per clus- 
ter, the individual cluster redshifts are uncertain to within 
700km/s (Dalton et al. 1994a). 



2.2 Simulations 

The cluster simulations are those generated by Croft & Efs- 
tathiou (1994a, b), using a particle-particle-particle-mesh N- 
body code (Efstathiou & Eastwood 1981, Efstathiou et al. 
1985), and are the same as those used by Dalton et al. 
(1994b). We consider three different models. First, the stan- 
dard CDM model which is scale-invariant and has fio = 1 
(hereafter T = fi/t = 0.5, CDM). Second, a scale-invariant, 
low-density CDM universe with fio = 0.2, h = 1 and a cos- 
mological constant A such that A/3Hq = 1 — fio (T = 0.2, 
CDM). The power spectrum for these two models is taken 
from Efstathiou, Bond & White (1992). Third, a scale- 
invariant model containing both a massive neutrino com- 
ponent (7-eV neutrinos) and CDM, such that Sltotai = 1-0, 
SIcdm = 0.6, Q v = 0.3, with h = 0.5 (Mixed Dark Matter, 
MDM) from Klypin et al. (1993). The models have been nor- 
malised to be consistent with the amplitude of fluctuations 
in the microwave background (Wright et al. 1994), so that 
the value of the linear theory rms mass fluctuations in 8 h -1 
Mpc spheres, as = 1.0 for the two CDM models and 0.67 for 
MDM. 10 realisations of each model were run, using different 
random phases, in boxes of side length 300 ft -1 Mpc, with 

6 

10 particles. We have also run 4 realisations of the T = 0.2, 
CDM with the slightly higher normalisation as = 1.3. Clus- 
ters of particles are located using a percolation algorithm, 
the cluster mass is defined within a metric radius and the 
richness limit of the cluster sample is chosen to match the 
observed cluster space density. The metric radius is cho- 



3 CLUSTER CORRELATIONS 
3.1 Method of estimation 

We apply a variant of the method used by Efstathiou et al. 
(1990), which accounts for selection effects automatically by 
estimating the mean density and its fluctuations in shells 
of mean constant redshift. We divide space in the redshift 
catalogue into a series of concentric shells, of width 2R, cen- 
tred on the observer, which are further subdivided into M s 
spherical cells of radius R. For each shell s we compute the 
mean number of particles N S (R) in a cell and the moments 
mj(R): 

m J = Y / (Nc-Ns) J , (1) 

c=l 

where N c is the number of galaxies in the c th cell. The cor- 
responding connected moments kj = N s £j, corrected for 
Poisson shot-noise (see section 3.4), are given by: 

k 2 = m 2 — N s 

ks = ms — 3wi2 + 2N s . 

ki = m± — 3m 2 — 6ms + llm 2 — 6N S . (2) 

(see e.g. Gaztanaga 1994). Thus, f 2 = k 2 / N s is the variance 
(also called a 2 ), f 3 = ks/N s is the skewness and f 4 = ki/N s 
the kurtosis. 

To estimate the mean value of (f 2 ) for the whole sample, 
we introduce the standard Gaussian likelihood function, so 
that the value of k 2 in each shell is weighted by the inverse 
of its variance, Var(k2) s : 

Var(A; 2 ) s = — , (3) 

where we have made the assumption that the cells are in- 
dependent. Thus the weighting compensates for the effects 
of fluctuations from both shot-noise (important for distant, 
diluted shells) or from having a small number of cells (impor- 
tant in the nearby shells). We have tried two prescriptions 
to estimate this variance. The first one uses the higher order 
moments mj in the shell to estimate kj in equation (3). In 
the second prescription we assume a Gaussian distribution 
to estimate kj from m 2 [as in Efstathiou et al. (1990)]. Both 
give similar results. 

For higher order moments we use the same weights for 
each shell as the ones used to estimate f 2 so that, in general: 
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Figure 1. The variance ^2^) estimated from counts-in-cells in 
a simulated mock catalogue of clusters (symbols with error-bars), 
compared with the values from clusters selected from the same 
realization of the simulation in a box (lines). The dashed line 
and filled squares correspond to the standard Q = 1 CDM model 
whereas solid lines and open squares correspond to low density 
CDM. 



M s 4 

= k g 4fe h ' (4) 

" M s 4 1 _1 

K = V Ns 

^ Var(k 2 ) s 

where { j = N^kj with kj evaluated for the s shell. 

One of the most valuable merits of this approach is that 
it also provides an estimate of the variance on this mean, 
Var(fj), as we have all the values in each shell and their 
relative weights: 

where M is the total number of shells, and the averages use 
the same weightings as equation (4). 

This is repeated for different values of R, the end prod- 
uct being the mean £j(R) and its variance. The cells in each 
shell must not overlap too much; this is important because 
otherwise fluctuations in M s are not taken into account in 
the weighting and nearby shells tend to dominate the statis- 
tics. 

3.2 Check of the method using simulated clusters. 

We have used the simulations described above to check our 
method. We first estimate the "true" values of £j(R) using 
counts in cells in redshift-space in the simulation box, in 
the standard way (see e.g. Baugh, Gaztanaga & Efstathiou 
1995). Next we transform the simulations to mimic the ge- 
ometry of the observational data, applying the APM survey 



Figure 2. The skewness £ 3 (i?) estimated from counts-in-cells in 
a simulated mock catalogue of clusters compared with the full 
simulation, as in Figure 1. 

mask and the selection function obtained by Dalton et al. 
(1994b). The above counts in cells method is used to esti- 
mate in the mock survey and compare them with results 
from direct estimation. 

In Figure 1 we show the comparison for mock T = 0.5 
and T = 0.2 CDM catalogues against the results from the 
original simulated box. In both cases we have included red- 
shift distortions. Figure 1 shows that there is good agree- 
ment, although it should be noted that a perfect match is 
not expected at large scales because the mock catalogue only 
samples a fraction of the total volume in the box. These re- 
sults are for just one simulation (not the average of 10 real- 
izations), showing that it is possible to recover the intrinsic 
clustering from the catalogues using the above method. 

Figure 2 shows the recovered skewness, £ 3 , from the 
T = 0.5 and T = 0.2 CDM catalogues compared with the 
results from the full cubical box. As expected the errors are 
much larger here. The errors in the T = 0.5 model are even 
larger and it is difficult to obtain significant results from 
just one catalogue given that the clustering amplitudes are 
smaller. 

3.3 Results for APM clusters and dark matter 
models 

Using the above method, we are able to measure correla- 
tions with a very large cell radius, up to R ~ 80 ft -1 Mpc. 
The smallest scales are limited by the shot-noise in the cell 
ocupancy to R ~ 5 ft -1 Mpc. The error bars come from equa- 
tion (5), which provides a consistent estimate of the variance 
of { j within the sample. 

To make an accurate estimate of { j in the simulations 
we use all clusters in each cubical box of side 300 ft -1 Mpc. 
There are ~ 1000 clusters in each simulation. In this way the 
uncertainties are typically smaller than the symbols we use 
to compare with the data. As mentioned above we have also 
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Figure 3. The variance ^2^) estimated from counts-in-cells in 
the APM clusters (symbols with with error-bars) compared to the 
values in different simulations. The continuous, short-dashed and 
long-dashed lines correspond repectively to the estimation in real 
space, redshift space and redshift space with Gaussian velocity 
errors. The lines show the standard Q = 1 CDM model (lower 
lines), the low density CDM variant (middle lines) and the MDM 
model (upper lines). 



made mock catalogues of the simulations, with the APM 
mask and selection function, and find a good agreement, 
within the errors, between the results from the mock cata- 
logue (using the shell method) and the direct results from 
the cube (see § 3.2 above). 

In Figures 3-5 we show the values of { 2 (-^-)> an d 
f 4 (i?), obtained from the APM clusters (symbols with error 
bars). The estimates for J = 4 become very uncertain, espe- 
cially at large scales. The "real space" results for { j in each 
model are shown as continuous lines in Figures 3-5. We es- 
timate { j for each model after applying redshift distortions 
using the peculiar velocities for the clusters in the simula- 
tions, short-dashed lines in figures. Finally, as the observa- 
tions are affected by errors in the redshift, we also show, as 
long-dashed lines, the estimates of in redshift space after 
adding errors with a Gaussian distribution of dispersion 700 
km/s to the simulated cluster redshifts. 

Redshift distortions produce slightly larger clustering 
amplitudes at large scales, because of coherent infall, and 
smaller amplitudes at small scales, because of random veloc- 
ities. The crossing of the real space and the redshift space 
curves indicates the scale at which both effects are compa- 
rable. We expect Gaussian errors of 700 km/s to reduce the 
redshift space results at small scales, ~ 700/_ffo, in agree- 
ment with the figures. 

In Figures 6 and 7 we show the hierarchical ratios Sj = 
?j/?2 with errors generated from the errors in fj. The 
mean values for the APM cluster are around S3 ~ 2 and 
S4 ~ 8. At the largest scales both the data and the models 
have large errors, nevertheless we show the results at these 
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Figure 4. The skewness £ 3 (i?) (symbols with error-bars) esti- 
mated from counts-in-cells in the APM clusters compared to the 
values in different simulated models, as in Figure 3. 
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Figure 5. The kurtosis £ 4 (i?) estimated from counts-in-cells in 
the APM clusters (symbols with error-bars) compared with the 
values in different simulated models, as in Figure 3. 



scales to illustrate that even when R ~ 50 ft -1 Mpc there is 
still a significant detection of S3 . 

Note that the simulations show that redshift-space dis- 
tortions produce only a small effect in Sj as the distortions 
in { j tend to compensate for the distortions in f 2 (Fry & 
Gaztanaga 1994). 

From Figure 3 it is apparent that f 2 in the T = 0.2 
model has too little power on scales R < 20 ft -1 Mpc. This is 
also evident in the skewness £ 3 , shown in Figure 4. Although 
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Figure 6. The hierarchical skewness 53 = f° r * ne APM 

clusters (symbols with error-bars) compared with values from dif- 
ferent cosmologicalmodels. The lower panel shows the low density 
CDM model, while the upper panel shows the MDM model. The 
continuous, short-dashed and long-dashed lines correspond repec- 
tively to the model estimation in real space, redshift space and 
redshift space with Gaussian velocity errors. 
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Figure 8. The variance C2 (^) m the ^ = 0.20 model normalized 
to <T8 = 1.0 (short-dashed line) and <T8 = 1.3 (long-dashed line) 
compared to ^2^) m the ^ = 0.15 model. 
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Figure 7. The hierarchical kurtosis S4 = C 4 /C 2 f° r the APM 
clusters (symbols with error-bars), compared with the values from 
different models, as in Figure 6. 
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Figure 9. The skewness ^(-R) in the V = 0.20 model normalized 
to <T8 = 1.0 (short-dashed line) and to <T8 = 1.3 (long-dashed 
line) compared to ^2^) m the ^ = 0.15 model (normalized so 
that (is = 1.0). 



we believe this effect is real, it is not very significant if we 
use 2 — a errorbars. 

While the amplitude of f 2 is not very sensitive to the 
normalization of the model, as, it does depend strongly on 
the value of T. To show this we have run one realization of 
a CDM T = 0.15 model (with Q = 0.2 and h = 0.75). This 
is shown in Figures 8 and 9, where we compare the results 
of the T = 0.15 and T = 0.20 models, with two normaliza- 
tions. Because we only have one realization of the T = 0.15 



model, these estimates correspond to just one set of random 
phases (identical for each model). With just one realization 
it is difficult to distinguish between these three models on 
scales larger than 20 ft -1 Mpc. The results at intermediate 
scales have sufficiently small errors to constrain the shape 
parameter T to have a value of ~ 0.15. This is in good 
agreement with results from the APM galaxy angular cor- 
relation function (Maddox et al. 1990a, Efstathiou, Maddox 
& Sutherland 1992). 
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3.4 Discreteness 



In both the observed cluster catalogues and the cluster sim- 
ulations the distribution of clusters is represented by a dis- 
crete number of points. The main effect of this discreteness is 
to introduce additional correlations (sometimes called 'shot- 
noise') which are important at scales R < d„, where d n is 
the mean interparticle separation. Typically, if we have a 
correlation length Ro, which is defined as £ 2 (i?o) = 1, with 
Ro > d„, as in the case of galaxies, then the intrinsic clus- 
tering dominates over the shot-noise. For clusters, Ro < d„, 
so the shot-noise dominates the measured clustering. 

This problem becomes evident when estimating £ j from 
moments of counts in cells because one has to decide explic- 
itly what (if any) discreteness correction should be applied. 
But the problem is equally important when estimating the 2- 
point correlation function, {2, directly. In the latter case one 
counts distances between pairs in the sample and compares 
them with similar counts in a random (Poisson) distribu- 
tion of points. Thus, the Poisson shot-noise correction has 
been widely used for clusters even when this is not always 
explicitly stated. 

One could choose not to correct the measured corre- 
lations for this discreteness as long as the same density of 
points is used in the comparison with models and other cat- 
alogues. This is also the case if one regards the cluster distri- 
bution as intrinsically discrete and richness dependent, i.e. 
characterised by its mean density and therefore by d„. 

On the other hand one could try to correct for discrete- 
ness and estimate the intrinsic correlations of an underlying 
density field. This estimate would be independent of the 
interparticle separation and would be more useful in com- 
parisons with models or other catalogues.^ The cluster dis- 
tribution can then be imagined as a continuous field which is 
traced by our particular point distribution. This continuous 
field is not the matter distribution, but a biased transforma- 
tion of it. As in the case of the galaxy distribution one can 
further assume that the location of the point drawn from the 
underlying fluctuations is the result of a (Poisson) random 
process (with probability proportional to local density). In 
this sense, the Poisson shot-noise model can be applied to 
correct for the discreteness. 

In our case, the Poisson correction is larger than unity 
at scales R < 1 9 ft -1 Mpc, as the mean cluster density, N , 
in spheres with R ~ 19 ft -1 Mpc is N = 1. This can be 
seen in Figure 10, where the uncorrected values of £ 2 for 
the T = 0.5 model (filled squares) and the T = 0.2 model 
(filled triangles) are compared with the values of £ 2 after the 
discreteness correction, with open squares (triangles) for the 
r = 0.5 (r = 0.2) CDM model. 

In the same Figure we have also plotted the Poisson 



shot-noise correction £ 2 



N oc R as a dashed line. 



Note that, for the T = 0.5 model, the shot-noise contribu- 
tion is larger than the intrinsic clustering contribution at all 
scales. Although we know that the intrinsic clustering prop- 
erties are different in these models, we can see in Figure 
10 how the uncorrected values of £ 2 approach the Poisson 



t One would then be able to disentangle, for example, the ef- 
fects of incompleteness (such as random sampling) and intrinsic 
richness dependence in a catalogue of clusters. 
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Figure 10. The uncorrected values of £ 2 for the T = 0.5 model 
(filled circles) and the T = 0.2 model (filled triangles) compared 
with the values of £ 2 after the discreteness (Poisson shot-noise) 
correction, with open circles (triangles) for the T = 0.2 (T = 0.5) 
CDM model. The Poisson shot-noise correction £ 2 = N oc R~ 3 
is shown as a dashed line. 



correction on small scales, where discreteness is important. 
Given that we have not used the Poisson model in deciding 
how to select clusters in the simulations this result indicates 
that Poisson shot-noise model is indeed a good approxima- 
tion to the discreteness effects. 



4 DISCUSSION AND CONCLUSIONS 

We find two common features in the cluster simulations of 
different dark matter models. First, redshift distortions are 
not very important. This is in spite of the fact that different 
models have quite different velocity fields (e.g. Croft & Ef- 
stathiou 1994b). Second, the values of S3 and S4 are quite 
similar in all models, i.e. in Figures 6 and 7. This is in spite 
of the fact that different models have quite different values 
of the variance, skewness or kurtosis (Figures 3-5). 

Both of these features can be understood in terms of 
biasing. Clusters are high peaks in the galaxy (mass) dis- 
tribution and they are therefore biased tracers of the mass 
(Kaiser 1984). On large scales one expects coherent stream- 
ing motions to enhance redshift space fluctuations (Kaiser 
1987), but the relative effect is suppressed (by a factor b, 
the ratio of the mass value to the cluster value of £ 2 ) for 
a biased distribution. On the other hand the amplitudes of 
Sj do not seem to be much affected by redshift distortions 
for either the mass distribution, according to perturbation 
theory (Bouchet et al. 1994), or for galaxies, according to 
observations (Fry & Gaztanaga 1994). This might explain 
why for £ j = Sj £ 2 , with J > 2, redshift distortions are 
also not very important. 

Why should Sj be similar in all models? Following Fry 
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& Gaztanaga (1993), a local bias relating galaxy or matter 
fluctuations 8 to cluster fluctuations 8 C : 

CO 

S c (x) = f[S(X)] = J2 b ^ k (X), (6) 

k = 

where b\ = b, produces amplitudes for the cluster distribu- 
tion 53 and 54 of: 

5 3 = b~ 1 (S 3 , m + 3c 2 ) 

5 4 = 6- 1 (5 4 , m + 12c 2 53, m +4c 3 + 12^) (7) 

where c 2 = 62 jb and C3 = 63 jb are the respective quadratic 
and cubic contributions from the biasing transformation in 
equation (6), and 5j >m are the hierarchical amplitudes in 
the underlying matter distribution. Here we have used the 
fact that both the galaxy distribution in the APM and the 
matter distribution in the simulations follow the hierarchi- 
cal relation ~ 5j{ 2 (see e.g. Gaztanaga 1994, Baugh, 
Gaztanaga & Efstathiou 1995). From the above expression 
we see that the underlying amplitudes 5j >m are typically 
suppressed by 1/6 when selecting high peaks. That the val- 
ues of 53 and 54 we find are quite similar in all models, 
even when the matter amplitudes S? Jim and 54 >m are signifi- 
cantly different (see Juszkiewicz, Bouchet & Colombi 1993, 
Gaztanaga & Baugh 1995) indicates that the biasing contri- 
bution in equation (7) is significant. For example, the typical 
variation between models of S? Jim for the mass distribution 
at large scales is A53 >m ~ 2. Everything being equal, the cor- 
responding variation in the cluster distribution, from equa- 
tion (7), is A53 < 0.6 for the typical value b ~ 3. Given the 
uncertainties in the measured values of 53 it is not surprising 
that this small variation is undetected. 

Recently, Plionis & Valdarnini (1995) and Cappi & 
Maurogordato (1995) have estimated 53 and 54 from an- 
gular and redshift positions of Abell and ACO clusters. The 
hierarchial ratios shown in Figures 6 and 7, 53 ~ 2 and 
54 ~ 8 are similar to their redshift-space results, but are 
slightly smaller than the ones found by Cappi & Mauro- 
gordato in the angular distribution, 53 ~ 3 and 54 ~ 15. 
It is not clear to us how significant this discrepancy is, as 
Cappi & Maurogordato (1995) use bootstrap errors, which 
are unreliable (see Baugh, Gaztanaga & Efstathiou 1995). 
On the other hand uncertainties from the deprojection of an- 
gular into 3-D statistics are not taken into account by Cappi 
& Maurogordato. This can introduce large errors given the 
uncertainties in the selection function of clusters. 

Plionis et al. (1995) have made predictions based on 
use of the Zel'dovich approximation to evolve cosmological 
density fields before selecting clusters. They rule out the 
T = 0.2, CDM model at a high significance level on the 
basis that it produces a cluster distribution with too high 
a variance and skewness on scales r > 20 /t -1 Mpc. In Ap- 
pendix A we have carried out tests on this method of simu- 
lating the cluster distribution and we find it to be unreliable. 
Problems arise from the inability of the Zel'dovich approx- 
imation to predict accurately the location of the compact 
objects which we are using to trace the density field. In the 
case of T = 0.2 CDM we have also run Zel'dovich reali- 
sations using the same parameters (particle number, grid 
spacing, box length and cluster richness) as Plionis et al. 
and have used the same method to evaluate the moments 
(from the cluster point distribution smoothed onto a grid 



with a Gaussian kernel, radius R). We find that the vari- 
ance and skewness are largely overestimated in comparison 
with the N-body results. In our test, the comparison is in 
real space and the discrepancy is ~ 50% for i?=20 and i?=30 
/t -1 Mpc. This suggests that although Plionis et al. use re- 
sults from a different observational catalogue of clusters (the 
combined Abell-ACO catalogue) to make their quantitative 
comparison, their conclusions are different and in our view 
less trustworthy than those which would be reached using 
full N-body simulations. 

In summary, we have shown that in the case of clus- 
ter samples discreteness has an important effect on clus- 
tering statistics at all scales. The Poisson shot-noise model 
has been widely used implicitly when estimating the 2-point 
correlation function for clusters. Here we have argued that 
the Poisson shot-noise correction should indeed be applied 
when estimating correlations {j in the cluster distribution. 
We reproduce previous conclusions by Dalton et al. (1994b) 
based on the 2-point correlation function, mainly that the 
standard Q = 1 CDM model is inconsistent with the second 
order statistics at all scales. As a new result we find that this 
model is inconsistent with either the third or fourth order 
statistics at all scales. A Mixed Dark Matter model gives 
good agreement, as does a low density CDM variant, both 
of these models being consistent with the { j observations 
for J = 2 — 4 at the 2a level. We have also shown that it is 
possible to measure 53 in the cluster distribution at scales 
as large as R ~ 50 ft -1 Mpc. The hierarchical correlations 
found in the cluster distribution are consequence of similar 
relations in the underlying matter field. In a forthcoming 
paper, we use the relations between the cluster values of 
Sj and the correponding values for the matter distribution 
[e.g. equation (7)] to investigate the biasing process in more 
detail (in preparation). 

Acknowledgements 

We would like to thank George Efstathiou for helpful 
comments on a previous version of this manuscript. EG ac- 
knowledges support of a Fellowship from the Commission 
of European Communities and from CSIC (Consejo Supe- 
rior de Investigaciones Cientificas) in Spain. RACC acknowl- 
edges receipt of a PPARC studentship. 

REFERENCES 

Abell, G.O. 1958, ApJ, Supp. 3, 211 

Abell, G.O., Corwin, H.G., Olowin, R.P. 1989, ApJ, Supp. 
70, 1 (ACO) 

Bahcall, N.A., Soneira, R.M., 1983, ApJ, 270, 20 
Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 

1986, ApJ, 304, 15 
Baugh, CM., Efstathiou, G., 1993, MNRAS, 265, 145 
Baugh, CM., Gaztafiaga, E. Efstathiou, C, 1995, MNRAS, 

in press 

Bernardeau, F., 1994a ApJ, 433, 1 
Bernardeau, F., 1994b A&A, 291, 697 
Bernardeau, F., Kofman 1994 ApJ, in press 
Bower, R.G., Coles, P., Frenk, C.S., White, S.D.M., 1993, 
ApJ, 405, 403 

Bouchet, F.R., Strauss, M.A., Davis, M., Fisher, K.B., Yahil, 
A., Huchra, J. P., 1993, ApJ, 417, 36 



8 Gaztanaga, Croft & Dalton 



Bouchet, F.R., Colombi, S., Hivon, E., Juszkiewicz, R., 1994, 

A&A, submitted 
Cappi A., Maurogordato, S. 1995, ApJ, 438, 507 
Coles, P., Melott, A.L., Shandarin, S.F., MNRAS, 260, 75 
Colombi, S., Bouchet, F.R., Schaeffer, R., 1993, Astron. As- 

troph., 281, 301 
Croft, R.A.C., Efstathiou, C, 1994a, MNRAS, 267, 390 
Croft, R.A.C., Efstathiou, C, 1994b, MNRAS, 268, L23 
Dalton, G.B., Efstathiou, C, Maddox, S.J., Sutherland, 

W.J., 1994a, MNRAS, 269, 151 
Dalton, G.B., Croft, R.A.C., Efstathiou, C, Sutherland, 

W.J., Maddox, S.J., Davis, M., 1994b, MNRAS, 271, 

L47 

Dalton, G.B., Efstathiou, C, Maddox, S.J., Sutherland, 

W.J., 1995, in preparation. 
Efstathiou, C, Bond, J.R., White, S.D.M, 1992, MNRAS, 

258, lp 

Efstathiou, C, Eastwood, J.W., 1981, MNRAS, 194, 503 

Efstathiou, C, Davis, M., Frenk, C.S., White, S.D.M., 1985, 
ApJ Supp.57, 241 

Efstathiou, C, Kaiser, N., Saunders, W., Lawrence, A., 
Rowan- Robinson, M., Ellis, R.S., Frenk, C.S., 1990, 
MNRAS, 247, lOp 

Efstathiou, C, Sutherland, W.J., Maddox, S.J., 1990, Na- 
ture 348, 705 

Efstathiou, C, Dalton, G.B., Sutherland, W.J., Maddox, 

S.J., 1992, MNRAS, 257, 125 
Frieman, J. A., Gaztanaga, E. 1994, ApJ, 425, 392 
Fry, J.N., Peebles, P.J.E. 1978, ApJ, 221, 19 
Fry, J.N. 1984, ApJ, 279, 499 
Fry, J.N., Gaztanaga, E., 1993, ApJ, 413, 447 
Fry, J.N., Gaztafiaga, E., 1994, ApJ, 425, 1 
Gaztanaga, E., 1992, ApJ, , 398, L17 
Gaztafiaga, E., 1994, MNRAS, , 268, 913 
Gaztanaga, E., Frieman, J. A., 1994, ApJ, 437, L13 
Gaztanaga, E., Baugh, CM., 1995, MNRAS, in press 
Goroff, M.H., Grinstein, B., Rey, S.J., Wise, M.B., 1986, 

ApJ, 311, 6 

Groth, E.J., Peebles, P.J.E., 1977, ApJ, 217, 385 
Hauser, M.G., Peebles, P.J.E., 1973, ApJ, 185, 757 
Hockney, R.W., Eastwood, J.W., 1981, Computer Simula- 
tions Using Particles, McGraw-Hill, New York. 
Jing, Y.P., Zhang 1989, ApJ, 342, 639 
Jing, Y.P., Valdarnini 1991, A&A, 250, 1 
Kaiser, N., 1984, ApJ, 284, L9 
Kaiser, N., 1987, MNRAS, 227, 1 

Klypin, A. A., Kopylov, A.I., 1983, Sov. Astron. Lett. 9, 41 
Klypin, A. A., Holtzmann, J., Primack, J., Regos, E., 1993 
ApJ, 416, 1 

Juszkiewicz, R., Bouchet, F.R., Colombi, S. 1993 ApJ, 412, 
L9 

Maddox, S.J., Efstathiou, G., Sutherland, W.J., Loveday, J., 

1990a, MNRAS, 242, 43p 
Maddox, S.J., Efstathiou, G., Sutherland, W.J., Loveday, J., 

1990b, MNRAS, 243, 692 
Maddox, S.J., Efstathiou, G., Sutherland, W.J., Loveday, J., 

1990c, MNRAS, 246, 433 
Mann, R.G., Heavens, A.F., Peacock, J. A., 1993, MNRAS, 

263, 798 

Meiksin, A., Szapudi, I. & Szalay, A.S. 1992, ApJ, 394, 87 
Nichol, R.C., Collins, C.A., Guzzo, L., Lumsden, S.L., 1992, 
MNRAS, 255, 21p 



Peebles, P.J.E., 1980, The Large Scale Structure of the Uni- 
verse: Princeton University Press 
Plionis, M., Valdarnini, R., 1995, MNRAS, in press 
Plionis, M., Borgani, S., Moscardini, L., Coles, P., 1995, 
ApJ, in press 

Romer, A.K., Collins, C.A., MacGillivray, H., Cruddace, 
R.G., Ebeling, H., Boringer, H., 1994, Nature, 372, 
75 

Sathyaprakash, B.S., Sahni, V., Munshi, D., Pogosyan, D., 
Melott, A.L., 1994, MNRAS, submitted 

Sharp, N.A., Bonometto, S.A., Lucchin, F. 1984, ApJ, 130, 
79 

Szapudi, I., Szalay, A.S. & Boschan, P. 1992, ApJ, 390, 350 
Toth, G., Hollosi, J., Szalay, A.S. 1989, ApJ, 344, 75 
Weinberg, D.H., Cole, S., 1992, MNRAS, 259, 652 
Wright, E.L., Smoot, G.F., Kogut, A., Hinshaw, G., Tenorio, 

L., Lineweaver, C, Bennett, C.L., Lubin, P.M., 1994, 

ApJ, 420, 1 
Zel'dovich, Ya., B., 1970, A&A, 5, 84 



Variance, Skewness & Kurtosis from galaxy clusters 9 



APPENDIX A: CLUSTERS IN THE 
ZEL'DOVICH APPROXIMATION 

In this appendix we are concerned with testing some aspects 
of the reliability of theoretical predictions of cluster cluster- 
ing. In particular we compare our N-body simulations with 
simulations that use the Zel'dovich (1970) approximation to 
evolve the density field. Such simulations have been used 
as a computationally inexpensive way of making predictions 
about clustering and higher order moments of clusters for 
a wide range of models (Plionis et al. 1995 and references 
therein). 

Al Alternative methods of simulating clusters. 

The two-point correlation function of simulated rich clus- 
ters evolves weakly with time (Croft & Efstathiou (1994a)). 
This is because much of the clustering measured is statis- 
tical, in the sense that the clusters were born clustered, as 
would happen if they formed at the locations of high peaks 
in the initial density field (Kaiser 1984). Less time consum- 
ing methods than the use of full N-body simulations might 
therefore be adequate to model the statistical clustering 
properties of clusters. Analytical approaches based on high 
peak theory (eg. Mann, Heavens & Peacock 1993) predict 
the same trends as N-body simulations - the cluster correla- 
tion amplitude increases weakly with cluster richness, for ex- 
ample. To make accurate predictions which can be compared 
quantitatively with observational samples, though, it is cer- 
tainly better to identify clusters with peaks in the evolved 
density field. 

This is the view that has been taken by Plionis et al. (1995) 
in their study of cluster correlations. The dynamical scheme 
which they use is the Truncated Zel'dovich Approximation 
(TZA), which has been tested by several authors (e.g. Coles, 
Melott & Shandarin 1993, Sathyaprakash et al. 1994) and 
found to give a surprisingly good indication of the behaviour 
of gravitating matter even into the non-linear regime. It 
should be noted that agreement found with N-body sim- 
ulations was only good on scales much larger than the typ- 
ical radius of galaxy clusters. We also note that it fails to 
reproduce in detail the amplitudes of the higher order cor- 
relations Sj of the mass at all scales (e.g. Juszkiewicz et al. 
1993). When the original Zel'dovich (1970) approximation 
is applied to a particle density field all particles move along 
their initial trajectories with constant velocity. When par- 
ticles trajectories cross one another, structure which would 
have formed under the action of gravity is erased. For this 
reason, in the improved version, the truncated ZA, the initial 
power spectrum is truncated (by smoothing with a gaussian 
filter) to remove the power at small wavelengths which is 
responsible for the bulk of orbit crossing. Using the TZA 
to evolve the density field rather than a full N-body code is 
computationally much cheaper, meaning that large numbers 
of models can be run and tested. This is potentially of bene- 
fit, as for example Plionis et al. (1995) have run simulations 
of 6 different universes, with which it would be interesting 
to compare our APM cluster statistics. 

It is therefore important to test the reliability of this sort 
of scheme, particularily as we have reached the stage where 
cluster samples are large enough that we can measure higher 
order statistics with relatively small errors. We would like to 



know whether clusters do indeed form at sites that can be 
predicted by weakly non-linear evolution of the density field. 
If this does occur, then it could give support to the idea that 
statistical clustering (or "high-peak biasing") is mainly re- 
sponsible for the distribution of clusters. On the other hand, 
if we are not able to find a reliable correspondence between 
TZA and N-body clusters the matter distribution may still 
be similar for both methods when smoothed on large scales, 
but in the TZA we will have been unable to predict cor- 
rectly the locations of clusters that we are using to trace 
that matter distribution. 

A2 Comparison of Truncated Zel'dovich 
Approximation and N-body clusters. 

We test which of these is in fact the case by evolving the mat- 
ter distribtion for the same initial random phases using the 
TZA and the P 3 M N-body code. In this example case we will 
use the T = 0.2 CDM models normalised so that <r& = 1.3. 
The input power spectrum of the TZA run was truncated 
by use of a Gaussian filter of radius r g = 6.3 /t -1 Mpc which 
Plionis et al. (1995) find is optimal for this particular model. 
The final particle distribution for a 15 /t -1 Mpc thick region 
is shown in Figure Alb(ii), and can be compared directly 
with the N-body particle distribution (Figure Ala(ii)). It 
can be seen that there are broad similarities between the 
two density fields. There are however fewer obvious clumps 
in the TZA density field. We are expecting this as small scale 
power has been removed from the power spectrum and also 
no mechanism exists within the dynamical scheme to form 
structure into virialised objects. We must therefore decide 
where such objects would have formed. Plionis et al. assume 
that this would occur at peaks in the TZA density field and 
we will follow them in defining clusters to be at the positions 
of peaks. However, examination of the TZA particle distri- 
bution shows that there are many smooth-looking regions 
which we would like to break up in order to reproduce the 
N-body results. 

To find clusters we grid the final particle density fields us- 
ing a TSC scheme (Hockney & Eastwood 1981) onto a 256 3 
mesh and then apply Gaussian smoothing. The positions of 
the 1000 highest peaks in the smoothed field are taken as 
the locations of the 1000 richest clusters. This technique does 
have a number of disadvantages, such as the fact that reso- 
lution is lost below the scale of the grid, and also due to the 
smoothing procedure. There is an additional anti-correlation 
effect, arising from the fact that peaks cannot lie in adjacent 
grid-cells. This is discussed further below when we calculate 
clustering statistics. To check the sensitivity to differences 
in cluster-finding techniques of the N-body results, we ap- 
ply the above procedure first to the N-body outputs. We try 
filters of radius r/ = 0.5, 1.0 and 2.0 /t -1 Mpc. As the inter- 
particle separation in cluster cores is small, we should not 
have to smooth much in order to remove any discreteness 
effects (in fact using the unsmoothed grid does give good 
results). In figure Ala(ii) we compare the positions of these 
N-body clusters (for r/ = 1.0 /t -1 Mpc) with those identified 
using percolation. It can be seen that there is an extremely 
good correspondence between the two cluster catalogues, as 
we would expect. The panel Ala(ii) shows the distribution 
of particles taken from the 100 /t -1 Mpc square region indi- 
cated in the left hand picture. 
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Figure 1. A 15 /i -1 Mpc thickness slice through simulations of T = 0.2 CDM universe with a% = 1.3 run using (a) a P 3 M N-body 
code and (b) the truncated Zel'dovich approximation. The solid dots in panels a(i) and b(i) show the postions of clusters identified in 
these simulations as high peaks after smoothing the density with a Gaussian filter of radius 1.0 /i -1 Mpc. The open circles show clusters 
identified as peaks with Gaussian smoothing of with filter radius 2.0 h~ 1 Mpc in panel b(i) and percolation groups in panel a(i). In panels 
a(ii) and b(ii) we zoom in on the 100 /i _1 Mpc square region enclosed by the solid lines in a(i) and b(i), showing the distribution of 
particles in the N-body and Zel'dovich simulations respectively. 



We apply this cluster-finding scheme to the TZA simula- 
tion. Because the nature of objects to be identified with 
clusters is not obvious from the particle plot, we try sev- 
eral different filter sizes r f = 1.0,2.0 and 3.0 /t _1 Mpc. We 
then look for peaks in the smoothed field. The results for 
Tf = 1.0 /t _1 Mpc and 2.0 /t _1 Mpc are shown in Fig. Alb(i), 
where we plot those of the 1000 highest peaks which fall in a 
slice of thickness 15 ft -1 Mpc. Again it is instructive to com- 
pare with the N-body results. We can see that neither value 
of if produces a good correspondence between TZA clus- 
ters and the N-body clusters. It also seems that the obvious 
caustics which the TZA produces are places where maxima 
are preferentially selected, particularily when the smoothing 
radius is small. A larger filter scale Tf than with the N-body 
results should be necessary in order to avoid any particle 



discreteness problems, as the peak regions do not have such 
a high particle density. We have run TZA simulations with 
200 3 particles in the same sized box in order to check on 
the effects of varying mass and spatial resolution. Our re- 
covery of similar results to the tests with 100 3 particles sug- 
gests that the mismatch between TZA and N-body clusters 
is caused mostly by the TZA dynamical scheme itself being 
inadequate. 



A3 Clustering statistics and Zel'dovich clusters. 

It is possible though that the two sets of objects can give 
statistically similar results, on the large (> 10 ft -1 Mpc) 
scales used in comparison with observations, and on which 
the TZA should be most useful because the density field is 




Figure Al. The variance £ 2 { r ) °f clusters identified from simulations of aT = 0.2 CDM universe with a% = 1.3 run using a P 3 M N-body 
code and the truncated Zel'dovich approximation. The error bars represent the rms error on the mean calculated from the scatter between 
4 realisations. In the N-body results we show £ 2 ( r ) f° r clusters selected from percolation groups (solid line) and from peaks in the density 
field (open symbols) smoothed with a Gaussian filter of radius 0.5 , 1.0 and 2.0 /i _1 Mpc. The filled symbols correspond to the different 
results obtained after using Gaussian filters of radius 1.0, 2.0, and 3.0 /i _1 Mpc to smooth the TZA evolved density field. Clusters were 
identified as the 1000 highest peaks in this field, giving a mean intercluster separation of 30 /i _1 Mpc. 



more linear. We have calculated the variance f 2 (r) f° r the 
TZA and N-body clusters, the results being shown in Fig 
A2. At this point it is worth mentioning that the method 
of selecting clusters from a grid (which is necessary in the 
TZA case) affects the correlations on small scales. The anti- 
correlation introduced by the fact that peaks must be a min- 
imum of 2 grid cells apart will make f 2 (r) artificially low on 
small scales. The Gaussian smoothing radius r/ used on the 
evolved field before peak selection will also depress £ 2 ( r )- 
This effect is visible in the plot for the N-body results: us- 
ing larger if results in deviations from the percolation curves 
out to larger r. The TZA clusters should also be affected by 
both problems, but it difficult to decide from the plot be- 
cause we do not know what properties a TZA sample with 
no gridding or smoothing would have. We will therefore con- 
centrate on scales > 10 /t -1 Mpc in the rest of our discussion. 
The N-body values of £ 2 ( r ) calculated after using percola- 
tion and peaks to find clusters are almost identical. The 
value of Tf which gives the worst agreement, 2.0 /t -1 Mpc, 
is probably too large for the compact objects present in the 
N-body output and in any case using a grid and filter is 
not our method of choice for selecting N-body clusters. For 
the TZA clusters, however, the correlation strengths depend 



strongly on the smoothing radius used to select clusters. It 
would appear that we are selecting a different population of 
objects as the smoothing radius is changed. With a small ra- 
dius there is marked excess clustering, as those peaks in the 
initial density field which would collapse to form large single 
clusters are instead broken up into many smaller objects. Us- 
ing if = 2.0 gives the closest results to the N-body clusters, 
although £ 2 ( r ) i s s ^ m marginally too low. It is important to 
realise that when comparing a range of Zel'dovich models 
realistically with observations, we should have an a priori 
reason for choosing the smoothing scale. The best scale may 
depend on the shape of P(k) or the amplitude of fluctua- 
tions. We have seen here that varying Tf between 1.0 and 3.0 
/t -1 Mpc gives enormously different results for £ 2 ( r )- Also, it 
apparent from the plots of cluster positions that even with 
the best radius we are not selecting the same objects as in 
the N-body output. This may imply that other properties 
of the cluster density field are not reliably predicted by this 
method. 
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A4 Discussion. 

As there is no clear visual correspondence between the N- 
body and TZA clusters, the exact nature of the objects we 
are selecting in the TZA density field is unclear. When we 
use a small r/ or the unsmoothed field, then it is possible 
that we are affected by particle discreteness, although our 
tests with higher resolution TZA simulations show that this 
is not a major problem here. This suggests that some of the 
objects found in the TZA output are artifacts of the limited 
dynamics, as are for example the shells of matter that have 
undergone obvious orbit crossing. By smoothing the initial 
power spectrum, as we are forced to do, and because the 
TZA does not work well on small scales, we are losing infor- 
mation which we need to locate where virialised objects will 
form. This is related to the problem that the efficacy of the 
TZA will depend on the overall amplitude of fluctuations. 
For higher normalisations, we will need a larger truncation 
scale on the power spectrum to suppress the increased orbit 
crossing. The objects that form will correspond even less to 
those in the fully non-linear calculation. This implies that 
the TZA could perform better in comparison with N-body 
simulations of the other models with different power spectra 
and lower as, which we have not tested here. On the other 
hand, as the galaxy clusters which we are using to trace the 
density field are by definition compact objects, and proba- 
bly virialised, our theoretical predictions should be able to 
predict the location of these objects. The Zel'dovich Approx- 
imation, whilst reproducing some other aspects of the large 
scale evolution of the density field, is unable to do this. 
In summary, we find that we can reliably select clumps of 
matter which we identify as galaxy clusters from N-body 
simulations, with different selection methods producing al- 
most identical results. This gives us confidence that our the- 
oretical calculations of { j will be robust. We have also tested 
a method for generating simulated cluster catalalogues us- 
ing the Zel'dovich approximation, and compared the clusters 
selected using this method with those in our own fully non- 
linear calculation. This leads us to conclude that the use of 
clusters from the Zel'dovich approximation is not really suit- 
able when engaging in a quantitative comparison of theories 
with observations of cluster clustering. 



